%% Series 8 - Exercise 1
% Table of Contents
% 1. Load data
% 2. Identification of the design
% 3. Determine which model could be fitted with the data
% 4. Statistical quality of the design
%         Dispersion matrix
%         The composite design as a reference
%         Variance Inflation factors
%  5. Half effects
%     5.2 Significant coefficients
%         Reference distribution
%         Normal plot
%         Stepwise fit
%         ANOVA ( Type II)
%          Summary of Model Comparison
%         Interpretation of Lack of Fit
% 6. Canonical analysis
%     6.1 Fix point determination
%     6.2 Axes determination
% 6.3 Visualization with slice
%%
% The objective of this exercise is to train the analysis of a design before and after the measurement. A by-product is the competence about plotting data, fitting models, producing ANOVA tables, performing a canonical analysis.
% the tasks to realize are:
% Load the data
% Identify if it is a classical design
% Check what type of model can be fitted
% Check the statistical quality of the design
% Infer the coefficients for a linear response model, a linear response model with interaction and a quadratic model. For each situation check the p-values and the LoF
% Perform a canonical analysis and get some insight of the response function
%
%%
%% 1. Load data
%
% The data of all the exercises are in the datafile "data_doe.xlsx". You can load them with the help of the specific tool "load data" in the "Home" menu.
% Data is put in a table because it allows the presence of get some metadata and also different types of data.

ipdata=readtable('data_doe.xlsx',...
    'sheet','Deposition',...
    'Range','A5:F28',...
    'ReadVariableNames',true,...
    'ReadRowNames',true);
disp(data)

% There is then 23 measurements for 4 factors, and one response.

Nexp=size(data,1);

%
%%
%% 2. Identification of the design
%
% The structure of the table and the title of the exercise give hint that the design is a Doehlert design. A quick geometrical analysis will confirm the fist impression.
% Let's start by extracting the essay matrix from the data table and normalize it :

Essay=data{:,1:4};
Essay_norm=rescale(Essay,-1,1,"InputMax",max(Essay),"InputMin",min(Essay));
disp(Essay_norm)

% The function  convhull()  will be used to determine the data points defining a convex space.  The a hull around the data points is draw  in 2D selecting the factors 2 by 2 :

for k=1:3
    for l=k+1:4
        figure
        % determining the border points
        Index=convhull(Essay_norm(:,k),Essay_norm(:,l));
        
        % plotting the profile of the border
        plot(Essay_norm(Index,k),Essay_norm(Index,l),'b-',"LineWidth",2)

        % plotting the data points
        hold on 
           plot( Essay_norm(:,k),Essay_norm(:,l),'ok',"MarkerFaceColor","red")
        hold off

        title (['plot ', num2str(k), '-', num2str(l)])
    end
end

% We may observe profils (projections on a surface) close to hexagons. We may also observe the number of levels of each factor: We are dealing with a Doehlert design for four factors. Matlab do not have a routine for generating this design. You can download the routine doehlert() from Moodle.
% A Doehlert design for 4 factors counts 21 runs. Here we have 23 runs

% Find unique rows and their indices
[~, uniqueIdx, groupIdx] = unique(Essay_norm, 'rows');

% Count occurrences of each row
counts = accumarray(groupIdx, 1);

% Find rows that appear more than once
duplicateRows = uniqueIdx(counts > 1);

% Display the duplicate rows
disp('Duplicate rows:');
disp(Essay_norm(duplicateRows, :));

disp('Number of duplications:');
disp(counts(counts > 1));

% The center experiment has been replicated 3 times.
%
%
%% 3. Determine which model could be fitted with the data
%
% A Doehlert matrix having been identified, the models that can be fitted are
% y=a_o+\sum a_i x_i+ \epsilon (linear response)
% y=a_o+\sum a_i x_i+ \sum_{i<j}  a_{ij} x_i x_j +\epsilon (linear response and interactions)
% y=a_o+\sum a_i x_i+ \sum_{i \leq j}  a_{ij} x_i x_j +\epsilon (quadratic)
% For the 3rd degree, there is definitivelly not a sufficient number of experiments as a Taylor polynome of 3rd degree having 4+4+12+15 =35 coefficients.
% Let's build a model matrix X for the quadratic model.
% For 4 factors, the quadratic model counts 15 degrees of freedom:
% one constant a_o
% 4 main effects a_i
% \frac{4!}{2!\ 2!}=6 interactions coefficients
% 4 quadratic coefficients

X_doehlert=x2fx(Essay_norm,'quadratic');

% the labels corresponding to the coefficients of the second degree model
coeflabels={'a_o','a_1','a_2','a_3','a_4',...
    'a_{12}','a_{13}','a_{14}','a_{23}','a_{24}','a_{34}',...
    'a_{11}','a_{22}','a_{33}','a_{44}'};

%
%%
%% 4. Statistical quality of the design
%%
%% Dispersion matrix
%
% The analysis of the quality of a design starts with the computation of the dispersion matrix.

% dispersion matrix of the doehlert design
D_doehlert=inv(X_doehlert'*X_doehlert);

%%
%% The composite design as a reference
% The composite design is recognised as a good design for response surface. So let's compare the performance of the two designs.

% Computing the composite matrix as a reference
E_composite=ccdesign(4,'center',3);
X_composite=x2fx(E_composite,'quadratic');
D_composite=inv(X_composite'*X_composite);

%
% The figure below compare the diagonal elements of the two designs which are directly related to the size of the confidence intervals.
% Doehlert: 21+2 runs
% Composite: 25+2 runs
% We observe that the quality of the Doehlert design is lower than the quality of the composite design. For the interaction and the quadratic coefficients the ratio is between 10 and 5 with only 4 runs of difference.
%

% A bar chart of the diagonal elements of the dispersion matrix
figure
bar([diag(D_doehlert),diag(D_composite)])
title('Diagonal elements of the dispersion matrix for quadratic model')
legend('Doehlert','Composite')
xticklabels(coeflabels)
xlabel('Coefficients')
legend("Position", [0.65,0.8,0.2,0.1])

% Compute the common color scale limits
minValue = min([D_doehlert(:); D_composite(:)]);
maxValue = max([D_doehlert(:); D_composite(:)]);

% heatmap of the dispersion matrix
figure
imagesc(D_doehlert)
colorbar
title('Dispersion matrix for the Doehlert design')


% Set color limits to be the same for both plots
clim([minValue maxValue]); 

% heatmap of the dispersion matrix
figure
imagesc(D_composite)
colorbar
title('Dispersion matrix for the composite design')

% Set color limits to be the same for both plots
clim([minValue maxValue]); 

% Let's now have a look to the whole dispersion matrix:
% We may observe a dispersion matrix standard for a second degree model with diagonal element of medium quality and covariances between quadratic coefficients and the constant.
% The Doehlert design is economical but its quality stay poor. It is its versatility that make it very useful.
% By comparison, the corresponding composite design with 27 experiments would have a dispersion matrix of a lot better quality as shown in the second image where covariance are only with the constant.
%
%% Variance Inflation factors
% The inflation factors give an indication of the consequences of colinearities. It can be observed a few  differences with the composithe design which is slightly better for the interaction coefficients and  quite equivalent for the quadratic terms.

% correlation matrix of the Doehlert design
C_doehlert=corrcov(D_doehlert); 

% VIF matrix of the Doehlert design
VIF_doehlert=diag(inv(C_doehlert)); 

% correlation matrix of the composite design
C_composite=corrcov(D_composite);

% VIF matrix of the composite design
VIF_composite=diag(inv(C_composite));

% comparison bar chart
figure
bar([VIF_doehlert(2:end),VIF_composite(2:end)])
title('Variance inflation factors)')
xticklabels({'a_1','a_2','a_3','a_4',...
    'a_{12}','a_{13}','a_{14}','a_{23}','a_{24}','a_{34}',...
    'a_{11}','a_{22}','a_{33}','a_{44}',})
legend('Doehlert','Composite')
legend("Position", [0.16,0.75,0.2,0.1])

%
% All these elements show that the the quality of the Doehlert design is poor in comparison of the composite. When ever possible the latter schould be prefered.
%
%%
%%  5. Half effects
% Let's compute now the half-effects for a quadratic model with the least square algorithm
% \alpha= (X^T X)^{-1} X^TZ Y
%  the residuals
% \epsilon=Y-\hat{Y}= Y-X \alpha
% and the experimental variance
% \sigma^2=\frac{\epsilon \cdot \epsilon}{N_{exp}-1}
% These calculations, for the normalised essay matrix, can be sumerize in the following bar chart:

% LSF estimates of the coefficients
coef=D_doehlert*X_doehlert'*data{:,5}; 

% Residuals
epsilon=data{:,5}-X_doehlert*coef;

% Experimental variance
s2=sumsqr(epsilon)/(Nexp-1);
s=sqrt(s2);

% Bar chart of the half-effects
figure
bar(coef(2:end)/coef(1),"green")
title('Half relative effects (Normalized) ','FontSize',14)
label={'a_1','a_2', 'a_3','a_4',...
    'a_{12}','a_{13}','a_{14}',...
    'a_{23}','a_{24}',...
    'a_{34}',...
    'a_{11}','a_{22}','a_{33}','a_{44}'};
set(gca,'XTickLabel',label)
xlabel('Coefficients')

% +/- sigma range
hold on
plot([0 15],[s/coef(1) s/coef(1)],'r--',"LineWidth",2)
plot([0 15],[-s/coef(1) -s/coef(1)],'r--',"LineWidth",2)
hold off

% legend
legend('half-effect','\pm \sigma')
legend("Position", [0.16,0.8,0.2,0.1])

%
%% 5.2 Significant coefficients
%
%% Reference distribution
% As there is sufficient degrees of freedom, an experimental variance has been computed and then if necessary a first screening can be done with the \pm \sigma interval as shown on the bar chart of the half-effects with the two dashed red lines.
%
%% Normal plot
% A normal plot can also be used for the screening of significant coefficients:

figure
h=normplot(coef(2:end));
title('Normal plot','FontWeight','bold','FontSize',16)
xlabel('Sorted coefficients','FontSize',14)
ylabel('norminv(p)- theoretical quantile','FontSize',14)


% draw again the sigmoid based on the coefficients not retained

% get the abcisses (coefficients)
x=get(h(1),'Xdata')'; 

% model matrix for the coefficients not retained
x=[ones(7,1),x(6:12)];

% get the ordinates (probabilities)
y=get(h(1),'Ydata')';
y=y(6:12);

% strait line regression
alpha=inv(x'*x)*x'*y; 

% delete the other strait line
delete(h(2))

% draw again  the red line
set(h(3),'Xdata',[-4,2],'Ydata', [1 -4;1 2]*alpha) 


[Essay,index]=sort(coef(2:end));
text(h(1).XData+.25, h(1).YData+.05, label(index),'FontSize',16)

% The following coefficients have been retained :
% a_o, a_2, a_3, a_{13}, a_{14}, a_{34},a_{22},a_{44}
% The model is then
% \hat{y} = a_o+ a_2 x_2 +  a_3 x_3 +  a_{13}x_1x_3+ a_{14}x_1x_4 + a_{34} x_3 x_4 + a_{22} x_2^2 + a_{44} x_4^2
%
%% Stepwise fit
% An alternative is the routine stepwiselm():

mdlstw=stepwiselm(Essay_norm,data{:,5},'quadratic','lower','constant','Penter',0.03,'PRemove',0.05)

% The obtained model is the same as long as the input criteria are standard. We may observe an anomaly in the routine that get rid of the terms  a_1 and a_4 although the high p-values élevées and the 'lower' parameter defined as 'constant'.
%% ANOVA ( Type II)
% To perform a precise ANOVA, a term matrix can be defined and then a type II ANOVA algorithm is applied

% define the spec matrix
spec=[0 0 0 0
    0 1 0 0
    0 0 1 0
    1 0 1 0
    1 0 0 1
    0 0 1 1
    0 2 0 0
    0 0 0 2];

disp(spec)
% Fit the model with specific terms
mdl=fitlm(Essay_norm,data{:,5},spec);

% compute the ANOVA table for each coefficient of the model
A_mdl=anova(mdl,'component',2);
disp(A_mdl)
plot(mdl)


% A lack of fit can be performed with the option 'summary' of the anova routine. A p-value of 96% indicates that there is no detectable lack of fit.

A_summary=anova(mdl,'summary');
disp(A_summary)

%
% Let's check what would have been with a simpler model without quadratic term

mdl_lin=fitlm(Essay_norm,data{:,5},'linear');
disp(mdl_lin)
A_lin=anova(mdl_lin,'summary')
mdl_int=fitlm(Essay_norm,data{:,5},'interactions');
disp(mdl_int)
plot(mdl_int)
A_int=anova(mdl_int,'summary');
disp(A_int)

%
% Let's summerize the situation

% build a table with 

Tmodel=categorical({'Linear';'Interaction';'Quadratic'});
%Tsignificant={'2,3';'2,3';'2,3,13,14,34,22,44'};
Tsignificant={[2,3];[2,3];[2,3,13,14,34,22,44]};
Tpval=[0.004;0.01;5e-9];
TloF=[0.17;0.14;0.96];

Tsummary=table(Tmodel,Tsignificant,Tpval,TloF,...
     'VariableNames',{'Model','Significant Terms','Model p-value','LoF'});
format shortG
disp(Tsummary)
format short

%
%%  Summary of Model Comparison
%
% Three models were evaluated to analyze experimental data collected using a Doehlert design with four factors: a linear model, an interaction model, and a quadratic model. Each model demonstrated strong statistical validity, as indicated by the high p-values for the lack of fit test, suggesting that none of the models failed to represent the observed data systematically.
% Linear Model: The linear model identified the main effects of factors 2 and 3 as significant, with p-values of 0.14% and 1%, respectively. The overall model p-value was 0.4%, indicating a good fit to the data, and the lack of fit p-value was 17%, confirming the adequacy of the linear model for describing the general trends in the data.
% Interaction Model: By including interactions, the second and third factors remained significant (p-values of 0.4% and 2%), but the interactions are not significant, with p-values above conventional thresholds.
% Quadratic Model: The quadratic model incorporated additional terms, including quadratic effects for factors 2 and 4 and significant interactions among several factor pairs. This model had a strikingly low p-value of 5 \times 10^{-9}, suggesting it captured the data's variability exceptionally well. The lack of fit p-value was 96%, indicating no significant evidence of systematic deviation between the model and observed data. This suggests the quadratic model not only fits the data but does so without overfitting.
%%
%% Interpretation of Lack of Fit
% The lack-of-fit test yielded high p-values across all models, indicating that none of the models systematically failed to represent the data. This consistency may result from the efficient experimental design, low noise in the data, or the ability of even simpler models to capture the primary structure of the response. However, while all models are statistically acceptable, the quadratic model is preferable as it accounts for more complex relationships, including significant curvature and interactions, providing a richer and more accurate representation of the underlying response surface.
% While the linear  model offer valid and interpretable fit, the quadratic model provides the most comprehensive description of the data. Its ability to capture significant main effects, interactions, and curvature makes it the most appropriate choice for understanding the response and making predictions.
%
%% 6. Canonical analysis
% The objective of the canonical analysis is to gain insight in the geometry of the quadratic model :
%  \hat{y}\approx51 -3 x_2 -2.2 x_3 -3.2 x_1 x_3 + 1.5 x_1x_4 + 1.8 x_3x_4 -2 x_2^2 -3.4 x_4^2
%
%% 6.1 Fix point determination
% The first step concists in deteremining the fix point of the function
%  \vec{x}_{s}= -\frac{1}{2} A^{-1} \vec{a} \quad \text{ et } \quad y_{s}= a_{o}+ \vec{x}_{s} \cdot \vec{a} +  \vec{x}_{s} ^T A \vec{x}_{s}

% get the coefficients of the fit
coef=mdl.Coefficients.Estimate;

% select the constant
ao=coef(1);

% select the linear terms
a= [0;coef(2);coef(3);0]; 

% build the curvature matrix
A= [0           0 coef(4)/2 coef(5)/2
    0           coef(7)  0  0
    coef(4)/2 0 0 coef(6)/2
    coef(5)/2 0 coef(6)/2 coef(8)];

% compute the fix point
xs= -inv(A)*a/2;

% compute the value of the model at the fix point
ys= ao + a'*xs + xs'*A * xs;

% So the coordinates of the fix point are within the domain :

% Place the values within a table
T_fp=table([xs;ys], ...
    'VariableNames',{'Fix point'},...
    'RowNames',{'x_1','x_2','x_3','x_4','y_s'});
format bank
disp(T_fp)
%format short

%
%% 6.2 Axes determination
% The geometry of the function is determined by analysing the eigen vectors and the eigen values

% compute the eigen values and eigen vectors
[V,lambda]=eig(A);

% place the values in a table
T_eig=array2table([V;diag(lambda)'],...
    'VariableNames',{'X1','X2','X3','X4'},...
    'RowNames',{'x1','x2','x3','x4','lambda'});

disp(T_eig)

% The analysis of the eigen vectors shows that
% A first eigen vector \tilde{X}_1 is in the direction of  x_4 with an eigen value \lambda_1\approx-4
% A second eigen vector  \tilde{X}_2 is in the direction of  x_2 with an eigen value \lambda_2 \approx-2
% Two eigen vectors  \tilde{X}_3  and \tilde{X}_4 are more or less parallel to the plane  x_1 x_3 with an eigen value \lambda_3\approx -1 and \lambda_4\approx 1.6
% The canonical  model is y=52- 4 \tilde{X}_1^2 - 2  \tilde{X}_2^2 -  \tilde{X}_3^2 + 1.6  \tilde{X}_4^2 that is an hyperboloïd in a 4D space
% Interaction are present only for axes x_1, x_3 and x_4 , we may then visualize the function in this 3D space x_1x_3x_4.
% Let's start by visualizing cuts with the function slice
%% 6.3 Visualization with slice
%

% Defining the points where the function must be evaluated 
[X1,X3,X4]=meshgrid([-1:.2:1]); %

% the values of X2 is fixed
X2=0.57; 

% Compute the function at the points of the network

Y=ao+ a(1)*X1 + a(2)*X2 + a(3)*X3 + a(4)*X4 +...
    A(1,3)*2*X1.*X3 +A(1,4)*2*X1.*X4+A(3,4)*2*X3.*X4+...
    A(2,2)*X2.^2 + A(4,4)*X4.^2;

% Planes cross the fix point

xslice=.78;
yslice=0.08;
zslice= [.16];

% Figure

figure
slice(X1,X3,X4,Y,xslice,yslice,zslice)
title('Value of the function within the experimental space')
xlabel('x1')
ylabel('x3')
zlabel('x4')
colorbar
map=colormap;
caxis('manual')
view([-24.70 40.40])

% The hyperbolique behaviour is easy to observe in the plane  x_3 x_4
% Keep in memory the scale for the color

Ncouleur=size(map,1);
[Cmin,Cmax]=caxis;

%
% Draw the iso-surfaces

figure
valeur= [44:2:54];
for k=1:6
    p=patch(isosurface(X1,X3,X4,Y,valeur(k)));
    isonormals(X1,X3,X4,Y,p)
    p.FaceColor=map(10+k*7,1:3);
    p.EdgeColor='none';
    hold on
end
hold off
daspect([1 1 1])
view(3);
axis tight
camlight
lighting gouraud
view([-48.30 8.40])
title('isosurface of y')
xlabel('x_1')
ylabel('x_3')
zlabel('x_4')

% We may now observe the complex behaviour of the response which is not straitforward visible looking at the function.
% Changing the value of x_2 (the fourth factor) between -1 and +1 we may observe the the behaviour stays the same with maximal values at the middle of the axes x_1=1, x_3=-1.
% The next step of the analysis should be with the subject matter for an scientific interpretation
% The point to keep in mind is the importance of the canonical analysis to understand and interpret the function of the second degree.